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Abstract: 

We  develop  a  parametric  statistical  approach  to  detecting  and  estimating  the  delay  and 
duration  parameters  of  a  ripple-fired  signal.  Such  signals  are  generated  by  mining  blasts 
with  roughly  equally  spaced  charges.  Detection  of  the  echo  structure  induced  by  ripple 
firing  can  serve  as  a  basis  for  discriminating  between  mining  blasts  and  nuclear  explosions 
or  earthquakes  which,  in  principle,  will  not  contain  such  echo  effects. 

The  model  assumes  a  noise  corrupted  signal  that  can  be  written  as  a  convolution  of 
source,  path  and  instrument  responses  with  a  pulse  sequence.  The  underlying  responses  are 
modeled  by  low  order  autoregressive  moving  average  (ARMA)  models  that  are  consistent 
with  conventional  deterministic  formulations  of  source  theory  and  instrument  response. 
The  pulse  sequence  is  modeled  as  a  seasonal  ARMA  process,  with  the  period  corresponding 
to  the  firing  delay  and  with  a  model  order  that  is  proportional  to  the  signal  duration.  We 
use  the  cepstrum  to  suggest  a  range  of  values  for  the  duration  and  delay  parameters  and 
then  search  all  possible  seasonal  ARMA  models  within  this  range.  The  final  model  chosen 
is  the  minimizer  of  the  corrected  Akaike  Information  Theory  Criterion  {AICc)- 

Limited  simulations  are  given  to  show  that  both  the  duration  and  delay  can  be  es¬ 
timated  effectively  with  the  seasonal  ARMA  search  when  the  delays  are  equal  and  that 
reasonable  estimators  are  available  when  the  delays  are  variable.  The  methodology  is 
applied  also  to  the  P  phase  of  a  single  Scandinavian  mining  explosion  where  we  obtain 
reasonable  estimators  for  delay  and  duration  under  the  assumption  of  ripple-firing. 

Key  Words:  Mining  explosions,  duration  and  delay,  echo  detection,  cepstra,  seasonal 
ARMA,  model  selection. 
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OBJECTIVE 


Regional  seismic  monitoring  and  discrimination  capabilities  that  are  desirable  imder  a  po¬ 
tential  Comprehensive  Test  Ban  Treaty  (CTBT)  can  be  improved  by  developing  algorithms 
and  new  procedures  for  distinguishing  between  earthquakes,  nuclear  explosions  and  mining 
explosions  of  various  kinds.  Much  effort  in  past  discrimination  studies  has  concentrated  on 
extracting  various  features  of  the  spectrum  that  are  characteristic  of  earthquakes,  nuclear 
explosions  or  mine  blasts. 

One  particular  spectral  feature  that  characterizes  some  mining  explosions  is  periodic 
modulation  of  the  spectrum  introduced  by  ripple-firing.  A  ripple-fired  event  involves  det¬ 
onation  of  a  number  of  mining  explosions  that  are  often  regularly  grouped  in  space  and 
time.  Such  explosions,  known  as  quarry  blasts,  have  low  magnitudes  that  may  be  close  to 
those  of  nuclear  explosions  that  one  might  monitor  under  the  CTBT.  Hence,  procedures 
for  detecting  ripple  firing  can  serve  as  bases  for  discriminating  between  mining  blasts  and 
nuclear  explosions  or  earthquakes.  A  number  of  authors  have  examined  various  aspects  of 
this  problem  and  have  proposed  techniques  for  analyzing  these  ripple-fired  seismic  signals 
(see  Baumgardt  and  Ziegler,  1988,  Hedlin  et  al,  1990,  Chapman  et  al  1992). 

The  objective  of  this  study  is  to  develop  robust  statistical  procedures  for  detecting 
and  characterizing  the  echo  structure  induced  by  ripple-firing.  This  includes  developing 
maximum  likelihood  in  both  the  time  and  frequency  domains  for  the  purpose  of  estimating 
the  parameters  corresponding  to  the  delay  and  duration  of  a  ripple-fired  signal  as  well  as 
the  spectrum  representing  source,  path  and  instrument  responses  that  will  be  convolved 
with  the  ripple-fired  signal.  We  use  modern  model  selection  techniques  based  on  Akaike’s 
Information  Theoretic  Criterion  (AIC)  to  determine  the  number  of  pulses  and  their  dura¬ 
tion. 


PRELIMINARY  RESEARCH  RESULTS 


We  have  developed  a  time  domain  approach  to  the  problem  of  detecting  and  identifying 
the  echo  parameters  associated  with  a  ripple-fired  event.  The  underlying  model  assumes 
that  a  sequence  of  n  signal  pulses  St  separated  by  d  will  produce  an  observed  process  of 
the  form 

n 

yt  =  Y^  ocjSt-jd  +  rit  (1) 

j=i 

where  rit  is  an  additive  noise.  The  observed  data  is,  therefore,  characterized  in  terms  of  a 
series  of  scaled  replicas  of  an  underlying  signal  that  last  for  a  total  duration  (n  -|-  l)d  and 
are  observed  at  delays  d,  2d, . . . ,  nd.  The  signal  St  is  an  underlying  process,  assumed  to 
satisfy  a  low-order  autoregressive  moving  average  model  of  the  form 


p 

k=l 


=  wt-  y^^dkWt-k- 

fc=i 
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(2) 


The  characteristics  of  the  signal  process  are  determined  by  the  p  autoregressive  parameters 
■  ■  ■  i4>p,  the  moving  average  parameters  ^2?  •  •  •  j  and  the  variance  of  the  white 
noise  process  Wt-  In  the  usual  Box- Jenkins  approach  (see  Shumway,  1988,  Chapter  3)  the 
model  determined  by  (2)  is  identified  as  ARMA{p,q).  It  is  convenient  to  rewrite  (1)  and 
(2)  using  the  backshift  operator  B,  defined  as  Bxt  =  xt-i,  and  we  obtain 


yt  -  otd(B)st  -b  nt 

(3) 

and 

4>(B)st  =  0(B)wt, 

(4) 

where 

n 

(5) 

i=i 

p 

(6) 

3=1 

and 

(7) 

i=i 

Now,  substituting  (4)  into  (3),  we  obtain 

yt  =  Xt  +  nt 

(8) 

where 

(f>{B)xt  =  d{B)ad{B)wt. 

(9) 

The  observed  process  yt  is  exhibited  in  (8)  as  the  sum  of  noise  and  a  multiplicative  ARM  A 
process  with  a  seasonal  moving  average  cic^(B)  of  order  n  and  season  d.  The  process  Xf 
is  identified  as  the  multiplicative  seasonal  ARM  A,  ARM  A{p,q)  x  (0,  n)^.  Under  this 
model  scenario,  the  process  of  detecting  a  ripple-fired  signal  is  reduced  to  looking  for  a 
multiplicative  seasonal  moving  average  component  in  an  ARMA  process,  where  the  order 
corresponds  to  the  number  of  pulses  and  the  season  corresponds  to  the  delay. 

The  problem  of  motivating  the  model  defined  by  (3)  and  (4)  can  be  approached  by 
appealing  to  Dargahi-Noubary  (1995),  who  shows  that  most  deterministic  source  models, 
including  that  of  Von  Seggem  and  Blandford  (1972)  can  be  generated  by  a  stochastic  model 
of  the  form  (4)  with  p  =  2  or  p  =  3  corresponding  to  the  and  uj^  models  respectively 
if  we  take  the  operator  as  (1  —  <pBy.  We  shall  not  specialize  the  operator  in  hopes  that 
the  more  general  form  (6)  will  help  to  adjust  for  instrument  response  and  path  effects. 
Furthermore,  these  effects  can  be  further  mitigated  by  adding  general  moving  average 
components  of  order  q. 

Although  the  model  suggested  by  (3)  and  (4)  is  not  exactly  ARMA(p,q)  x  (0,  n),; 
because  of  the  additive  noise  term,  it  is  well  known  that  ARMA  signal  plus  noise  models 
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are  equivalent  to  pure  ARM  A  models  with  additional  moving  average  components  so  we 
are  justified  in  restricting  candidates  to  models  of  the  multiplicative  form.  What  is  required 
is  estimating  the  ARMA  parameters  <^1, . . . ,  (j)p,  6i, . . .  ,6g,  as  well  as  the  amplitudes  of 
the  reflections  ni, . . . ,  This  can  be  done  for  a  fixed  p,  q,  n  and  d  using  nonlinear  least 
squares  on  the  residuals,  minimizing 


(10) 


where  Wf  can  be  obtained  by  solving  (9)  successively,  using  yt  in  place  of  Xf  There  are 
generally,  an  imlimited  number  of  values  of  p.,q,ri  and  d  possible  so  we  need  a  method 
for  deciding  which  values  are  most  reasonable.  A  common  model  selection  criterion  is  the 
corrected  form  of  Akaike’s  information  criterion  (see  Hurvich  and  Tsai,  1989), 


AICc  =  log  + 


T  +  pA  q  +  n 
T  —  p  —  q  —  n  —  2' 


(11) 


where  we  choose  the  model  that  minimizes  AICc-  For  the  large  sample  used  here,  this 
measure  gives  the  same  result  as  the  conventional  AIC . 

In  order  to  illustrate  the  above  methodology,  we  consider  analyzing  two  contrived 
events  and  a  single  mining  explosion  from  Scandinavia  shown  in  Figures  1,  2  and  3.  The 
contrived  events  were  simulated  by  generating  an  ARM A{2,0)  process  of  the  form 


St  -  st-i  +  .6st_2  =  wt  (12) 

for  the  signal  with  cr^  =  1.  This  is  Equation  (2)  with  p  =  2,  <^i  =  1,  ^2  =  -“-6.  The  process 
was  added  and  delayed  using 

yt  =  atXt  +  nt,  (13) 

where 

+  St-8  +  +  St-24  +  •St-32  (14) 

which  is  just  a  delayed  signal  with  n  =  4,  d  =  8,  i.e.  a  signal  with  four  pulses  separated  by 
8  points,  leading  to  a  duration  of  32  points.  To  make  the  output  resemble  an  explosion, 
the  ARMA  process  was  modulated  by  a  function  of  the  form  at  =  Bi  exp{— 02^}  with 
logdi  =  —7,^2  =  .01.  The  second  contrived  event,  shown  in  Figure  2  adds  five  pulses  at 
irregular  delays  and  noise  rit  with  variance  =  (-01)^.  Figure  3  shows  a  P-phase  from 
a  presumed  mining  explosion  recorded  at  FINESS  in  Scandinavia  and  given  by  Blandford 
(1993)  in  his  Table  1. 

To  analyze  these  three  examples,  we  search  for  the  best  ARMA(p,q)  x  (0,  n)^  model 
using  AICC  and  the  following  general  guidelines.  The  operators  (6)  and  (7)  that  best 
emulate  the  convolution  of  source,  path  and  instrument  response,  i.e.  p  and  q  in  the  model, 
are  restricted  to  low  orders.  In  this  case,  we  took  p  =  2, 3,  g  =  0, 1,  corresponding  to  values 
suggested  by  conventional  source  models  and  response  functions.  These  low  orders  all 
do  a  respectable  job  of  emulating  the  observed  spectra.  In  order  to  limit  somewhat  the 
search  over  the  seasonal  components  n  and  d,  the  spectra,  log  spectra  and  cepstra  were 
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Figure  1.  Analysis  of  a  contrived  signal  composed  by  modulating  an  ARMA(2,0)  process  and  delaying 
using  4  pulses  separated  by  8  points.  Frequency  is  in  cycles/point.  Delay  is  in  points. 
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Figure  2.  Analysis  of  a  contrived  signal  composed  by  modulating  an  ARMA(2,0)  process  and  delaying 
using  5  irregularly  spaced  pulses  at  intervals  8,  8,  10,  8  and  10  points  with  additive  noise.  Frequency  is  in 
cycles/point.  Delay  is  in  points. 
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Figure  3.  Analysis  of  a  P-phase  from  a  presumed  mining  explosion  recorded  at  FINESS  (FlAl)  in  Scandi¬ 
navia  (latitude  59.476,  longitude  24.1442)  dated  12/10/91  with  local  magnitude  2.59  (see  Blandford,  1993, 
Table  1).  Frequency  is  in  cycles  per  point  at  40  sec~^  so  that  the  folding  frequency  corresponds  to  20  Hz. 
Delay  is  in  points  or  .025  seconds. 
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computed.  The  cepstrum,  in  particular,  is  helpful  in  suggesting  values  for  n  and  d.  For 
example  see  Bogert  et  al  (1962)  who  suggested  the  cepstrum  or  Baumgardt  and  Ziegler 
(1988)  who  showed,  by  expansion,  the  periodicities  in  the  spectrum  expected  when  the 
spectrum  of  Xf  =  ad(B)st  is  computed.  Shumway  and  McQuarrie  (1994)  show  that  when 
ai  ~  a2  =  ■  ■  •  =  ajii  the  spectrum  of  xt  is 

sin^(7ruj(n  +  l)d)  ^ 
sm^{Trujd) 

where  Ps(<^)  is  the  spectrum  of  the  signal.  The  cepstrum,  i.e.  logP2;(a;),  then  will  be 
the  sum  of  the  logarithm  of  the  signal  spectrum  and  two  other  components  with  periods 
proportional  to  d,  the  delay  and  nd,  the  duration. 

To  begin,  consider  analyzing  the  first  contrived  series  shown  in  Figure  1  which  we  know 
to  be  approximately  ARMA{2, 0)  x  (0, 4)8  by  inspecting  the  generating  equations  (12)-(14). 
This  means  that  there  are  4  pulses  of  length  8  points  in  the  train  leading  to  a  duration  of 
5  X  8  =  40  points.  Both  the  spectrum  and  the  log  spectrum  show  periodicities  and  we  can 
see  the  two  distinct  components  in  the  log  spectrum.  Computing  the  cepstrum  shows  main 
peaks  at  8  and  40  corresponding  to  the  delay  and  duration.  To  check  the  seasonal  ARM  A 
search  procedure,  consider  searching  the  models  ARMA{2  —  3, 0  —  1)  x  (0, 3  —  8)5-12  which 
indicates  looking  at2<p<3,0<5<l,3<n<8, 5<d<12,  i.e.  low  order  ARM  A 
models  with  seasonal  components  between  5  and  12  and  orders  between  3  and  8.  Table  1 
shows  the  results  obtained  from  the  best  seven  models  under  the  column  Figure  1.  We  see 
that  the  minimum  corrected  AIC  is  for  the  correct  model. 

The  best  model  has  the  form  of  (12)-(14)  with  the  estimated  parameters  substituted, 
so  that 

St  -  .97st_i  +  .60sf_2  =  wt 

becomes  the  signal  model  and 


) 

(15) 


Xt  —  St  A  .99st_8  +  .87S(_i6  +  .80st_24  +  •71st_32 

gives  the  estimated  reflection  pattern.  Note  that  the  coefficients  are  not  all  near  imity 
because  of  the  exponential  modulation  of  the  Xt  in  the  observed  data  yt-  Table  1  also 
gives  the  next  best  six  models  and  we  note  that  all  but  one  has  the  correct  delay  but 
that  the  number  of  pulses  and  hence,  the  duration  changes.  The  spectrum  implied  by  this 
ARMA{2,0)  X  (0,4)3  model  is  shown  in  Figure  1  and  we  note  that  it  agrees  quite  well 
with  the  narrow  band  spectrum  (1  %  bandwidth). 
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Table  1.  Model  Summaries:  ARMA(p,q)  x  (0,  n)^ 


Figure 

P 

n 

d 

AICc 

1 

2 

0 

4 

8 

-7.2546 

1 

2 

0 

5 

8 

-7.2543 

1 

2 

0 

6 

8 

-7.2533 

1 

2 

0 

8 

8 

-7.2523 

1 

2 

0 

7 

8 

-7.2513 

1 

2 

0 

3 

8 

-6.8295 

1 

2 

0 

6 

6 

-6.2294 

2 

2 

0 

7 

8 

-5.7271 

2 

2 

0 

6 

8 

-5.7105 

2 

2 

0 

5 

8 

-5.7028 

2 

2 

0 

8 

6 

-5.6993 

2 

2 

0 

7 

6 

-5.6886 

2 

2 

0 

6 

9 

-6.6858 

2 

2 

0 

8 

9 

-5.6848 

3 

2 

1 

6 

6 

14.3649 

3 

2 

1 

4 

6 

14.3706 

3 

2 

1 

5 

6 

14.3714 

3 

2 

1 

3 

6 

14.3773 

3 

2 

1 

3 

12 

14.3830 

3 

2 

1 

4 

12 

14.3848 

3 

2 

1 

5 

12 

14.3865 

When  some  irregularity  in  the  pulse  intervals  is  introduced  as  in  Figure  2,  the  cep- 
strum  gives  quite  ambiguous  results.  From  Table  1,  note  that  the  seasonal  ARM  A  search 
indicates  a  model  with  7  pulses  separated  by  8  when  the  correct  model  has  5  pulses;  3 
are  separated  by  8  points  and  2  are  separated  by  10  points.  The  dmation  predicted  is 
(7  -f  1)  X  8  =  64  points  as  compared  with  the  known  duration  of  44  points.  The  ARM  A 
spectrum  shown  in  Figure  2  is  a  reasonable  version  of  the  narrow  band  smoothed  peri- 
odogram.  The  second  best  model  shows  6  pulses  separated  by  8  points  and  an  estimated 
diuration  of  56  points  which  seems  somewhat  closer  to  the  correct  model  as  is  the  third 
best  model  with  5  pulses  separated  by  8  points.  This  is  a  case  where  the  cepstrum  gives 
no  insight  into  the  structure  but  the  seasonal  ARM  A  search  finds  a  reasonable  model. 
Again,  note  that  the  spacing  is  easier  to  estimate  than  the  duration  because  the  number 
of  piilses  at  a  given  duration  all  have  about  the  same  value  for  AICq  and  other  durations 
have  generally  larger  values. 

In  Figure  3  is  shown  the  P-phase  from  a  single  mining  explosion  along  with  the  spectra, 
log  spectra,  cepstrum  and  ARM  A  spectrum  derived  from  the  best  seasonal  ARMA(2, 1)  x 
(0, 6)6  model.  The  sampling  rate  is  40  points  per  second,  leading  to  a  folding  frequency  of 
20  Hz  corresponding  to  .5  cycles  per  point  as  shown  on  the  graphs  in  Figure  3.  Table  1 
shows  that  the  four  best  models  all  have  a  duration  of  6  points  or  6/40  =  .15  seconds.  A 
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delay  of  150  milliseconds  (msec.)  may  be  somewhat  long  for  most  mining  explosions  (see 
Baumgardt  and  Ziegler,  1988,  Hedlin  et  al,  1990  or  Chapman  et  al,  1992),  but  it  is  clear 
that  this  sampling  rate  may  not  be  high  enough  to  pick  up  in  the  range  25-50  msec  and 
we  may  be  looking  at  aliases.  Hopefully,  there  will  be  data  available  at  higher  sampling 
rates  that  can  unequivocally  sort  out  reflections  in  the  desired  ranges. 

RECOMMENDATIONS  AND  FUTURE  PLANS 

A  high  priority  for  this  study  is  obtaining  more  multiple-phase  (P  and  S  arrivals)  explosion 
and  earthquake  data  sampled  at  a  higher  rate,  say  at  100  points  per  second.  We  would 
also  like  to  obtain  high  quality  signals  from  an  array  of  elements  so  as  to  be  able  to  take 
advantage  of  replication,  as  is  done  in  stacking,  practiced  by  both  Baumgardt  and  Ziegler 
(1988)  and  Hedlin  et  al  (1990).  For  this  to  work,  the  increased  difficulties  inherent  in 
replicated  ARM  A  searching  may  dictate  a  more  semi-parametric  likelihood  approach  in 
the  frequency  domain. 

There,  we  would  replace  the  ARM  A  model  for  the  spectrum  by  a  more  non-parametric 
version,  employing  the  general  form 

P,(u:)  =  |a4e-="“)pP,(u.)  +  (16) 

where  the  signal  spectrum  Ps{uj)  can  be  estimated  non-parametrically  using  the  Whittle 
likelihood  and  assuming  that  broad-band  smoothing  will  eliminate  the  duration  and  spac¬ 
ing  ripples.  We  may  also  be  able  to  assume  that  the  signal-to-noise  ratio  is  constant  over 
frequency  and  write 

P,i^)  =  P,(...)(ia,(e-^’"")P  +  iY  (17) 

where  =  Ps{u;) / Pn{uj)  is  the  signal  to  noise  ratio.  If  Ps{uj)  is  reasonably  constant  over  a 
broad  band,  it  can  be  replaced  by  the  smoothed  spectrum  over  that  interval  in  the  Whittle 
likelihood.  Such  likelihoods  can  be  treated  for  arrays  as  in  Der  et  al  (1992). 

In  the  frequency  domain,  we  may  also  consider  the  possible  of  taking  into  account 
the  amplitude  modulation  function  at  by  considering  the  Whittle  likehhood  of  the  time 
varying  spectrum  or  sonogram.  Dahlhaus  (1995)  has  shown  that  local  averaging  over  time 
of  the  Whittle  likelihood  can  be  used  to  estimate  parameters  of  the  modulating  function. 

A  final  note  relates  to  the  fact  that  the  spacing  values  d  will  not  be  the  same  and,  in 
fact,  may  be  quite  variable.  In  this  case,  we  must  replace  jd  in  (1)  by  dj  so  that 

n 

yt  =  Y^OijSt-dj  +nt,  (18) 

7=1 

where  di,  (i2,  •  •  • ,  form  a  sequence  of  unequal  delay  times.  Estimating  these  as  pa¬ 
rameters  seems  intractable  because  of  the  high  dimensionality  of  the  resulting  parameter 
space.  An  option  that  expect  to  try  is  letting  the  delay  times  be  random  and  possibly 
uniformly  distributed  on  the  integers  if  one  makes  no  specific  assumptions  about  the  time 
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delays.  One  may  also  assume  that  the  random  time  delays  are  clustered  about  the  points 
d,  2d,  3d, . . . ,  nd  in  some  specified  fashion.  The  integrated  likelihood  function  can  then  be 
computed  for  any  preset  combination  of  parameters  using  Monte-Carlo  methods. 
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